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Abstract. To study epitaxial thin-film growth, a new model is introduced and 
extensive kinetic Monte Carlo simulations performed for a wide range of fluxes 
and temperatures. Varying the deposition conditions, a rich growth diagram is 
found. The model also reproduces several known regimes and in the limit of low 
particle mobility a new regime is defined. Finally, a relation is postulated between 
the temperatures of the kinetic and thermal roughening transitions. 
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1. Introduction 

The development of experimental techniques to probe and visualize at the atomic 
scale has challenged studies on film morphology and the identification of the relevant 
controlling parameters. Several experimental [TJ [H [3J H] and theoretical [5J [5] 
studies have been performed to identify the underlying relevant mechanisms and 
processes at the atomistic scale. Herein, we introduce a model for multilayer growth 
to study the influence of deposition conditions, namely, temperature and deposition 
flux. 

Growth is an intrinsically nonequilibrium problem where each deposition event 
perturbs the system and keeps it out of equilibrium. It competes with thermally 
activated relaxation events that lead to equilibrium. Particles can arrive and stick 
to the substrate at a certain flux or hop between different basins, in the energy 
landscape, through thermally activated processes. Two different characteristic times 
are of relevance: the typical time between consecutive deposition events, which is 
related to the incoming flux of particles, and the hierarchy of relaxation times, related 
to temperature. The ratio between the first and the latter measures the departure 
from equilibrium [9]. 

During homoepitaxial growth, where deposited and substrate particles are the 
same, at temperatures below the roughness transition, the layer-by-layer growth 
regime is expected with nucleation of subsequent layers occurring after the previous 
one had time to complete. However, if the time between consecutive deposition events 
is short, on average, as compared to the typical timescale for the rate determining 
processes leading to relaxation, then thermal equilibrium is disrupted. Growth takes 
place away from equilibrium since thermal processes leading to relaxation remain 
unaltered. For example, one can expect island formation on top of existing islands, 
therefore leading to the nucleation of new layers before existing layers had a chance 
to become complete. 

We are primarily interested in the cooperative behavior that is established at long 
timescales and large length scales. We, therefore, depart from techniques that provide 
a more refined evolution at much shorter timescales and/or length scales [TU1 [TT1 [T2] . 
To this end we adopt a kinetic Monte Carlo approach, where the list of possible 
processes is identified a priori EH E3 EH S3- In this context, the properties of 
the model reflect the, a priori, chosen set of processes. Since this set of processes is 
known, it is easier to gain insight on the relevant mechanisms. 

There has been a substantial amount of effort focused on the study of 
submonolayer deposition, island formation and growth, and island-island coalescence 
[5J [T51 dni [201 1211 H2]> while for multilayer growth several questions remain open 
(23 [231 [251 Uni HZ] • Typically, either a detailed approach is considered, which attempts 
to make no assumptions, or as few as possible [H [28l ESI EOl [3_TJ [32l E3j, or simple 
models are developed, where mechanisms are restricted to the most relevant ones 
[MllMllMllSZllMllMliniST]. In this work, we adopt the latter approach, which 
allows us to better understand the relevance of each process. To rely on as few as 
possible parameters, we introduce a model where the activation barrier of each process 
is computed based on the one for the terrace diffusion and the coordination number 
in the initial and final states. This accounts for well known phenomena like, e.g., the 
corner rounding barrier. We acknowledge that the neighbor counting approach has 
known limitations, since pair-wise interactions are considered, thus neglecting many 
body effects. Despite the shortcomings, these models have successfully explained 
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relevant phenomena by addressing length and time scales of the order of the ones 
considered in this work (H US] 33] . 

Exploring the two-parameter space of temperature and deposition flux we not 
only reproduce several growth regimes previously reported, but we also identify a new 
one, namely, at low particle mobility. In this new regime, the typical analysis based 
on surface properties is not able to characterize it. Besides, the transition between 
the various identified regimes is also studied. 

The manuscript is organized in the following way. The model is introduced in the 
next section and simulation results introduced in section [3j along with their discussion. 
Concluding remarks are drawn in section [4] 

2. Model and definitions 

To study epitaxial growth we propose a model where particles are randomly deposited 
on the substrate, while deposited particles can then hop in the energy landscape 
between adjacent basins through thermally activated processes. The set of possible 
processes on the surface is given by the lattice sites surrounding the initial and final 
neighboring sites for a total of 2 10 configurations along each of the four possible 
directions, namely, up, down, left, and right. The activation barrier for each of the 
possible configurations is simply taken as the difference between the coordination 
number in the initial and final positions |32j . 

Ei = E m - (n/ - n ) J pp , (1) 

where E m is the activation barrier for terrace diffusion, J pp the energy of the particle- 
particle interaction, and the coordination number before (after) the hop is uq (n/). 
In spite of its simplicity, this strategy has been able to grasp relevant phenomena like, 
e.g., the high vacancy mobility observed in several experiments [44], and the relevance 
of perimeter mass transport for submonolayer island evolution on Ag(100) substrates 
[19] . The rate of process i is obtained from the Arrhenius equation, 

n = Vi e- E ^ T , (2) 

where is the attempt frequency. For simplicity the same attempt frequency is 
assumed equal for all processes, so Vi = v [45 . 

A key process in this system is the corner rounding |461 147) as illustrated in 
figure [T] This process is usually included by considering a direct jump to the step 
below, which requires an additional parameter for the activation barrier [29, 46 , 48, 49 . 
In the present model, the rules naturally account for such barrier, with a two steps 
process. 

The ballistic rule for deposition implies that deposited particles stick to the first 
occupied (nearest) neighbor, so particles can overhang [50] [51] [52] ■ However, the 
ballistic deposition model produces too many vacancies. Other models (see for example 
Ref. [381 [53]), such as solid-on-solid (SOS) and downward tunneling mechanism, 
have been considered for deposition 05], but they are not able to reproduce the 
experimentally observed concentration of vacancies. 

For epitaxial growth, particle detachment from the substrate can be neglected [54j . 
To achieve such a constraint, processes which break connectivity to the substrate are 
forbidden. To identify these situations we consider the Hoshen-Kopclman algorithm 
[55] , This rule prevents clusters from detaching from the substrate, a mechanism that 
can be neglected for a large region of the diagram of flux and temperature, and it its 
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a) 




Figure 1. a) Scheme of the corner rounding process and of the b) energy 
landscape. 



usefulness is solely restricted to the prevention of unphysical processes common on 
(KMC) lattice models. 

We introduce cj> as the ratio between the binding and thermal energies, 

<t> = lhr- ^ 

k B l 

The attempt frequency for the terrace diffusion process Vd is defined as 

_ Em 

Vd = ve k B T . (4) 

Now, if we also define the ratio, Re = E m /J pp , i.e., between the activation energy 
for terrace diffusion, E m , and the pair-wise interaction between particles, J ppi we can 
rewrite equation Q as 

v d = ve- R *+. (5) 

According to equations ^ and (|5j), the rate of a certain process i can be calculated 

by, 

n = v d e An ^. (6) 

We have casted the set of parameters of the model in terms of two, related to energy 
ratios, namely, Re and <f>. The former measures the activation barrier of the terrace 
diffusion in units of the particle-particle interaction while the latter provides a measure 
of mobility. The last parameter, the flux, F, defined as the number of impinging 
atoms arriving, per unit of surface (a line in the present case) and per unit of time, in 
monolayers per second (ML/s). 

In the next section we identify and analyze different growth regimes obtained 
with the present model and discuss their transitions. 



3. Results and discussion 



To reproduce growth regimes, we performed extensive kinetic Monte Carlo simulations 
of the (l+l)-dimensional system. The time increment is computed as At = 
— lnC/Si^ij where £ is an uniformly distributed variable in the range (0, 1] and 
the sum is over all processes accessible to each of the particles, thus giving the total 



Transitions between epitaxial growth regimes: A (l+l)-dimensional KMC study 5 



rate of the system. Simulations were carried out with periodic boundary conditions 
in the horizontal direction. Unless otherwise stated, simulations were performed on a 
lattice of 10 3 sites and results were averaged over 10 3 samples. 

We now start with the analysis of the effect of the flux of impinging particles 
in the morphology of the system. We compute the time evolution of the roughness 
defined as, 

W(t) = ^([h(x,t)-(h(t))} 2 ), (7) 

where h(x, t) is the height of the column x at time t and (h(t)) is the average height 
at time t. 

In figure [2^a) the evolution of the roughness is plotted for different values of the 
flux, with = 3 and Re = 5 kept fixed. After few deposited layers, different behaviors 
are observed and limiting regimes can be identified like, e.g., ballistic deposition at 
high fluxes and layer-by-layer at low deposition fluxes. For a characterization of the 
different growth regimes at various conditions the most commonly used measurement 
is the growth exponent /3, given by, 

W{t)~t . (8) 

In figure [2jb), we plot the value of /3 as a function of 1/F, for Re — 5 for different 
values of <fr. At fixed temperature (constant </>), the inverse of the flux is a measurement 
of the mobility in the system. At low mobility, the ballistic deposition (BD) regime 
is recovered, characterized by shadowing of lower layers due to the overhanging of 
particles. As the mobility increases, particles relaxation is promoted and the roughness 
increases; we denote this regime as ballistic deposition with local relaxation (BDLR). 
As the density of vacancies vanish, a plateau in the roughness exponent as a function 
of the flux is observed, since small islands are formed, at the level of the first layer, 
and interlayer diffusion favors relaxation |56j (at lower values of <fi the plateau is 
suppressed). At lower fluxes, the size of islands increases, which makes interlayer 
diffusion less likely to occur. Consequently, island nucleation atop of existing islands 
leads to stratified growth and an increase in the roughness exponent. This growth 
mode, where mounds are formed, corresponds to the kinetic rough (KR) regime. 
For even lower values of flux (higher mobility), the system has enough time to 
relax between deposition events and a layer-by-layer growth (LBL) is observed [51j . 
Snapshots of various regimes are shown in figure [3] 

In the following subsections, these different transitions and their growth regimes 
are discussed. First, the transition to the thermal rough regime is addressed, then 
the one from BD to BDLR, and finally the transition from BDLR to kinetic rough 
growth. This section ends with the characterization of the transition from kinetic- 
rough to layer-by-layer growth. 



3.1. Thermal roughening transition 

The roughening transition, or thermal roughening transition, as it is also known in the 
literature [8], intrinsically occurs under equilibrium conditions, though the remaining 
transitions referred to below occur away from equilibrium. To study this transition 
two different approaches are considered. We first consider its influence on the growth 
exponent j3 (figure EJa)) and we systematically study the behavior of an initially 
smooth surface at different temperatures (figure El b)). 
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Figure 2. a) Roughness evolution at <p = 3 as a function of time defined in units 
of deposited layers, for F = {5, 10 3 , 10 7 , 10 10 }ML/s. At higher fluxes the growth 
regime is rough and at lower fluxes the roughness oscillates representing a layer- 
by-layer growth, b) The roughness exponent /3 as a function of X/F for different 
values of <j>. From left to right <j> = {l(red), 2(green), 3(blue), 4(violet), 5(cyan)}, 
measurements performed after 10 2 deposited monolayers. 




abed 



Figure 3. Snapshots of the different growth regimes a) BD (<j> = 3; F = 1 X 10 15 ), 
b) BDLR O = 3;F = 4 X 10 6 ), c) kinetic rough (<f> = 3; F = 1 X 10 2 ), and d) 
layer-by-layer (<j> = 3; F = 1), for a system of fOO lattice sites and 100 deposited 
layers. Each color corresponds to 10 consecutively deposited layers. 



The transition to the layer-by-layer regime is observed for a wide range of <j> values, 
and occurs when the growth exponent becomes null. However, in figure |4]ja), below a 
certain value of 4>, the growth exponent does not converge to zero, i.e., the transition 
to a layer-by-layer growth does not occur. In fact, the roughness exponent /3 diverges 
at values of <fr < 2.0, which can be interpreted as a critical <p c . To systematically study 
this equilibrium transition we consider an initially flat surface with 10 layers and let it 
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Figure 4. a) Growth exponent evolution as a function of l/F for various <f> values 
near tj> c . From bottom to top we have <f> = {2, 1.8, 1.4, 1.2, 1.1}. b) Roughness as 
a function of the temperature for lattice sizes, from bottom to top, of 100, 200, 
400, 1000 and 2000 on a previously deposited film of 10 ML. 



evolve, for different values of 0, to a constant roughness, which required 6 x 10 4 to 10 s 
time units. In figure Qb), the roughness as a function of 4> for different system sizes 
is shown. The plot shows a logarithmic increase of the roughness below the transition 
point c/> c , as expected [S]- 

3.2. BD to BDLR 

At high fluxes or low temperatures, two growth regimes are present, BD and BDLR 
(see figure [3] (a) and (b), respectively), which differ by the presence of a partial 
restructuring during growth in the BDLR regime. However, this mobility is very 
low and particles mainly hop to nearby sites of higher coordination, i.e., sites with 
processes of higher activation barriers according to equation ([2]). 

Despite similar behavior of the roughness for both regimes, the snapshots in 
figure [3] reveal different bulk structures. Therefore, the typical analysis based on the 
roughness exponent [57] does not provide a way to differentiate between both regimes. 
Since the average coordination is higher in the BDLR regime, in figure[5|a) we analyze 
the behavior of the density of dangling bonds (bonds to perimeter sites) given by, 

where iVf, is the number of dangling bonds, L is the linear size of the system, and JVj 
is the number of deposited layers. 

In figure [5ja) a constant density of dangling bonds is observed for higher fluxes 
which decreases at lower fluxes. The point where the density of dangling bonds starts 
to decrease, depends on the mobility of the particles. In the thermodynamic limit 
one expects the restructuring, at finite values of the flux, to relentlessly lead the bulk 
of the deposit to a more compact, layer-by-layer morphology. However, our interest 
focuses on the roughness and the number of dangling bonds at the interface, which 
represent quantities that depend on the kinetic evolution of the growing interface. In a 
sense, the impact of mobility provides a measure of deviation from the extreme case of 
BD to BDLR. Consequently, the characterization of the transition is of a more subtle 
nature, based on the number of dangling bonds, while the growth exponent remains 
unaltered. For computational efficiency, we instead count the dangling bonds even in 
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Figure 5. Density of dangling bonds (a), pi,, and its variance, ipi, (b) as a 
function of l/F for different values of <j> after 100 deposited layers. From left 
to right, <f> = {1,2,3,4,5}. The insets are at <j> = 3 and for system sizes of 
L={400, 1600, 3200}. 



the bulk, since simulations are always carried over a finite number of (time) steps. At 
sufficiently low mobility one does not expect substantial restructuring of the bulk, but 
accounting for these bonds provides much better statistics. An interesting behavior is 
found for the variance of the density of dangling bonds ifb in figure [5^b) , given by 

W = ((p b -( Pb )) 2 ). (10) 

A maximum is observed, which can be taken to identify the transition from BD to 
BDLR. A non- vanishing maximum with the system size (inset of figure [5](b)) is a sign 
of a discontinuous transition in the thermodynamic limit |58j . Besides, the peak for 
the second moment shows a sharper discontinuity at the transition as particle mobility 
decreases. 

3.3. BDLR to kinetic rough 

As the flux lowers, particle diffusion becomes more relevant at a given temperature 
and particles placed on overhanging positions hop to more stable (higher coordination) 
positions, in part due to a larger relaxation time between consecutive deposition 
events. This behavior reduces the shadowing effect characteristic of the BD and BDLR 
regimes and, consequently, an increase in the roughness is observed (see figure [2](b)) . 
Since the number of overhanging particles is reduced, the density of vacancies in the 
film also decreases. We measure this density, given by, 

* - 1 - (id 

where Ni is the number of deposited layers plotted as a function of l/F in figure [6]j a). 
A transition on p v is found, where it has a value of 1/2, characteristic of the BD 
regime, and vanishes in the kinetic rough regime. As previously stated, the density of 
vacancies in BDLR does not significantly differ from that of BD, due to shadowing. 
Figure |6jb) shows the variance of the density of vacancies 

p v = ((p v -( Pv )) 2 ), (12) 

where (...) represents an average over samples. A maximum is observed, at the 
transition from BDLR to kinetic rough growth. However, the maximum is expected 
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Figure 6. The density of vacancies (a), p v , and its variance (b) as a function 
of 1/F for different values of <f> after 100 deposited layers. From left to right, 
<j> = {1,2,3,4,5}. The insets are at cj> = 3 and for different system sizes of 
L={400, 1600, 3200}. 



ln(F) 





Figure 7. Growth mode diagram with the kinetic roughening transition between 
layer-by-layer and kinetic rough growth, a) The flux as a function of <j> for values, 
from bottom to top, of Re = 10, Re = 5, and Re = 1. b) The ratio R as a 
function of cj> where triangles represent Re = 10, circles represent Re = 5, and 
squares represent Re = 1. 



to vanish in the thermodynamic limit, a sign of a continuous transition (inset of 
figure [6](b) ) |58j . A snapshot of the kinetic rough regime is shown in figure [3^c) . 

3.4- Kinetic rough to layer-by-layer 

Another relevant transition reproduced by this model is the kinetic roughening 
one between layer-by-layer and kinetic rough growth. The transition point can be 
computed by extrapolating in figure [2jh) the flux F, when the growth exponent, f3, 
goes to zero. A snapshot of the layer-by- layer regime is in figure [3jc) and|3^d). 

In figure[7^a), the transition points from kinetic rough to layer- by-layer growth can 
be observed. In agreement with previous work [511 152| , the transition has an Arrhenius 
dependency of the flux. For different values of Re (figure |7[a)), different slopes of 
the Arrhenius exponential function, representing the kinetic roughening transition, 
are observed. A ratio R is defined to relate the flux of particles F and the terrace 
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Figure 8. Final growth mode diagram, where TR and KR are the thermal rough 
regime and kinetic rough regime respectively. 



diffusion rate Ud, 

F 

R=—. (13) 

The plot of ln(i?) as a function of </>, figure^b), shows no influence from the various Re 
values, which represents, within the limits of the present study, a material independent 
property. Unlike the previous cases, for this transition a finite-size study has not been 
considered since the layer-by-layer regime represents a transient regime, which, at long 
timescales or large system sizes, crosses over to the kinetic rough regime [55] . 

A relation can now be established between the thermal roughening temperature 
T tr and the one from the kinetic roughening T^ r . This can be inferred from figure [7^b) , 
and the thermal roughening transition relation <j> tr = J pp /T tr as 

^ = ^exp(-^), (14) 

where A and a are independent of the temperature and flux, so that they solely depend 
on Re- This equation fits the plots of figure [7]Ja), where the Arrhenius dependency 
of the flux over the temperature is observed. To reduce the dependency on Re in 



equation (14), the ratio R (equation (13)) is used, leading to the relation 



R = Bexp[-c^ ), (15) 



T, 



At 



where B and c are constants. In figure [7](b) are the same lines of figure ^ a) with the 
proposed relation. 



3.5. Growth mode diagram 

The results are summarized in figure [HJ The diagram represents, for different pairs of <fi 
and F, the corresponding growth regime. BD/BDLR and BDLR/KR transition lines 
were picked from figures [SJb) and|6^b), respectively. For the KR/LBL and thermal 
roughening transition results from sections |3.4| and |3.1| were considered. The diagram 
of figure [8] represents the observed growth modes and the transition lines between 
them. 
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4. Concluding remarks 

We presented a KMC model to study the influence of both the flux and temperature 
on the epitaxial growth. We simulate an (1+1) -dimensional lattice mode with nearest- 
neighbor interactions that includes deposition and thermally activated processes. 
Despite the simplicity of the model, known regimes are reproduced like, e.g., ballistic 
deposition, kinetic rough, layer-by-layer, and thermal rough growth. Additionally, 
a new regime between BD and KR is identified, the ballistic deposition with local 
relaxation. 

We studied the transitions between the various regimes with emphasis on the 
BD/BDLR and BDLR/KR ones, where the traditional surface analysis has to be 
complemented by other properties like the number of dangling bonds or the number 
of vacancies. The transition KR/LBL was also studied and the Arrhenius relation 
between temperature and flux confirmed. Furthermore, we showed a relation between 
the thermal and kinetic roughening transitions. 

A more detailed study of the scaling properties is left open as well as the 
straightforward extension to a study in 2+1 dimensions. 
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